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The article proposes a bi-objective optimization approach for 
layup design of laminates. The optimization method 
combines the Particle Swarm Optimization (PSO) heuristics 
and Simulated Annealing (SA) optimization method. The 
minimum weight optimization is subjected to design 
constraints such as_ strength, stiffness, layup blending 
continuity, and several manufacturing design rules, which 
are combined as a single function and included within the bi- 
objective formulation. Several composite materials design 
problems are included to show the capabilities and 
usefulness of the proposed method. The optimization 
analysis has also been connected to the finite element 
analysis to solve the problem of composite plate optimization 
with blending constraints. The plate is divided into some 
regions, and the blending constraints are imposed globally by 
using the concepts of the greater-than-or-equal-to blending to 
achieve continuity of laminate layups across the regions. The 
results generally showed that the proposed method led to 
excellent results, representing a promising approach for the 
design of laminated composite materials. 
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1. Introduction 


Composite materials are part of the advanced materials with growing applications in many 
engineering fields. These materials have specific functions, and their directional properties 
enable to tailor them for a specific application. Complex mechanical behavior associated with a 
large number of variables causes the structural design of composite materials much more 
difficult and laborious than the structures made of isotropic (metallic) materials. These 
characteristics have encouraged the use of mathematical optimization methods in determination 
of the optimum shape and layup stacking sequence of the structure in a more sophisticated way 
to attain the maximum material efficiency and performance. Efficiency and performance of 
composite materials are in close relationship with design methodologies to define lamina 
thickness and orientation, resin and fiber materials, manufacturing, etc. to achieve reliability and 
resonable strength-to-weight and stiffness-to-weight ratios [1]. 


Research on the subject of composite materials optimization has been reported extensively. 
Generally, research in the field of optimum composite design focuses on the three areas of 
modeling, analysis, and optimization method, each of which is reviewed briefly. 


e Modeling 


The first step in the design optimization of laminated composites is the choice of material 
behavior modeling. Undoubtedly, the simplest way to optimize composites is when material 
behavior is expressed by closed-form mathematical relationships. If a closed-form solution 
exists, cost function computation and the sensitivity analysis will be managed efficiently. Since 
closed-form solutions are only available in simple cases and they generally involve many 
assumptions, this approach is normally conducted in the design of laminates at the material level. 
Employing the Classical Lamination Theory (CLT), effective stiffness parameters, stresses, 
critical buckling loads, and some other physical characteristics of laminates can be expressed by 
closed-form equations. In material level optimization analysis, the thickness of lamina, fiber 
orientation, the volume percentage of fibers, and other parameters in lamina level are commonly 
considered as design variables. For the sake of simplicity and computational efficiency, much of 
the research in the field of composite optimization has focused on material-level design. 
Optimum weight design of laminates with given stiffness properties [2], maximizing buckling 
load-carrying capacity of laminated composite [3], optimum design of laminates against fatigue 
loads [4] are only a few examples from literature, in which closed-form solutions are used for 
optimization purposes. For more complicated models such as reinforced plates, plates with 
cutouts, structures with 3D fiber placement, curved shells and 3D structures numerical methods 
such as finite element method is currently used [5]. 


e Analysis 


To date, several theories have been developed to determine kinematic behavior and stress 
distributions in composite laminates. Robbins & Reddy (1993) [6] categorized the analysis 
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methods based on the displacements or stress distributions through the laminate thickness to two 
classes: Equivalent Single Layer (ESL) and Layer-wise (LW) theories. The main difference of 
these methods is on the assumption of shear stress distribution in the thickness direction and 
transverse deformation of the laminate. 


ESL-based theories, reduce a 3D problem to a 2D one by defining the displacement as a linear 
combination of pre-defined functions and mid-plane displacement components. ESL is further 
classified into CLT [7], First Order Shear Deformation (FSDT) [8], and Higher-order Shear 
Deformation (HSDT) theories [9]. 


Classical lamination theory ignores shear deformations through the thickness and involves only 
three variables. However, in FSDT, shear deformation in the thickness direction is estimated by a 
linear function and therefore, it involves five variables. The main drawback of FSDT is that the 
equilibrium equations are not satisfied in the upper and lower surfaces of the plate. To 
compensate for the difference between the linear assumed stress and the real stress distributions, 
a corrective coefficient is often necessary. This coefficient is a function of the stack sequence and 
cannot be defined uniquely. In contrast, HSDT uses non-linear functions such as cubic, 
trigonometric, or hyperbolic functions to describe shear deformation in the thickness direction of 
the laminate. HSDT satisfy equilibrium equations but contains more variables, and thus is more 
computationally expensive [10]. Although ESL may be adequate for analysis of most practical 
composite laminates design problems, they typically fail to describe accurately the 3D stress 
field at the ply level. This deficiency is primarily related to the assumption of continuous 
transverse strain components across the interface of the lamina, which causes discontinuity of the 
transverse stresses. 


In LW theories, each layer is treated separately, and specific displacement field expansions or 
constraints are applied within each layer to yield a sort of continuity across the layer interfaces 
and a more accurate description of the complete stress state [11]. For instance, in a research 
work, each layer of a laminate was considered as an independent plate, and the compatibility of 
displacements at layer interfaces have been enforced by using the Legendre polynomials [12]. 


In the optimization process, different layups are examined iteratively, and hence, the calculation 
of the shear stress correction factor is crucial and costly. In this respect, for thin laminates, CLT 
is the easiest choice. However, for thick laminates, FSDT or other high order theories with an 
approximate shear stress correction factor, or finite element analysis should be used, even though 
it enforces much computational cost to the analysis. 


e Optimization methods 


Mathematical optimization methods, both deterministic and evolutionary stochastic search 
approaches, are widely used in the design of composite structures [13]. In the most design cases, 
the objective is the minimization of the total thickness or weight of the structure subjected to 
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constraints such as stiffness, strength, natural frequencies, critical buckling loads, and 
manufacturing requirements. 


Deterministic methods usually require objective function to be continuous and differentiable. 
Moreover, these methods can readily converge to a local optimum points. However, the 
application of these methods in real large design problems involving discrete and/or integer 
variables is limited. On the other hand, stochastic random search optimization methods can 
handle combinations of continuous and non-continuous functions, and can approach to the global 
optimum solution. However, evolutionary search methods generally necessitate a large number 
of objective function evaluation for thorough exploration of the feasible region. Besides, an 
algorithm is always needed to prevent convergence to local optimum points [14]. 


For design optimization of composite laminates, stochastic search methods have been addressed 
broadly in the literature. Genetic Algorithm (GA) [15,16], Simulated Annealing (SA) [2,6], Ant 
Colony [17], Particle Swarm Optimization (PSO) [18,19] and hybrid methods such as PSO-SA 
[20] and PSO-GA [21] are among the methods that are referred to extensively. 


In this article, a bi-objective optimization procedure based on the hybrid PSO-SA method is 
proposed for the design of laminated composites. The adopted hybrid PSO-SA optimization 
method is inspired from the method previously developed by Javidrad, & Nazari (2017) [22] and 
implemented to design of laminates by using penalty functions [20]. Stiffness, strength, layup 
blending continuity, and several design and manufacturing rules were incorporated into the 
formulation. Both CLT solutions and finite element modeling were used to compute the stiffness 
and strength properties of the laminate. Three design problems, including A) stiffness parameter 
maximization, weight optimization of a laminate under in-plane forces and moments, and weight 
optimization of a plate under bending loads including layup blending constraint were presented. 


2. Analysis of composite laminates 


Consider a laminated composite plate as shown in Fig. 1(a), with reference plane z=0 located 
at the mid-plane. Based on the Kirchhoff’s hypothesis, the displacement components are 


u=u pecs 
0 oz 
jay (1) 
Oy 
W=W, 


where u, v and w are small displacements along the x, y and z directions, respectively. 
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(b) 


Fig. 1. (a@) Geometry of a laminate with global (xyz) and material (1-2) coordinate systems. (6) Definition 
of force and moment resultants. 


The forces and moments per unit length of the cross-section can be found by integrating stress 
components in the thickness direction. Eqs. 2 describes the in-plane stress resultant forces [N,, 
Ny, Nxy]', and resultant bending moments [M,, My, Mxy]'. 


N, n om 
N, = » J ° Oy dz (2a) 
L Ny — xy 
[ M. n on 
M, |=>[. | 0, [za (2b) 
M.,, = xy 


where 6x, Gy and Oxy are two normal and shear stresses developed in lamina, respectively, and n 
represents the number of layers in the laminate. By implementing constitutive relations in global 
coordinates [23], the relationships between the resultant forces and resultant moments and the 
global mid-plane strains (€) and curvatures («) can be derived as below: 


[V]=[4][e]+[4][«] 
[M]=[8][e]+[P][«] (3) 


where A, B and D are stiffness matrices. The components of these matrices are calculated as: 
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Ay =>, QO) (Xa) 1b 7=12,6 (3a) 
k=l 
14, = be 

By = 5 Lj) (Zi — Zi) Lp = 12.6 (3b) 
k=l 
1S yes 

Dy = 5 2 (Qy)n Zi ~ 2-1) i, j =1,2,6 (3c) 
k=l 


where (O,) , denote components of the reduced stiffness matrix of the Ath lamina in the global 


coordinate system. If the same materials is used for each lamina in a laminate, the stiffness 
matrices can also be expressed in invariant forms [23]. 


2.1. Stress-based failure analysis 


In the design of laminated composite materials, a comprehensive strength-based failure criterion 
is required. The failure criterion must be able to predict the failure of a laminate for any layup 
configuration. However, due to the interaction between the layers, accurate failure prediction of a 
laminated composite materials is rather complicated [24]. 


To date, several strength-based failure criteria have been suggested. Limit-based theories 
(maximum stress or strain) [25], polynomial-based theories [26,27], strain energy-based theories 
(e.g. Tsai-Hill failure criterion)[28], and direct mode determining theories [29,30] are among 
those that extensively used in practice. Rohwer (2015) [31] presented a review of the failure 
criteria for fiber composite materials. 


In this study, failure is determined by the “first ply failure” (FPF) approach following the 
maximum stress (MS) and Tsai-Wu (TW) ply-level failure criteria. In many references (e.g. [32]) 
certain combinations of failure theories have been suggested for actual failure analysis to limit 
shortcomings of each failure criterion. 


A micromechanical study [33] showed that for a lamina with fiber orientation angles between 0° 
and 30° the TW criterion is accurate enough, while for fiber orientation angles between 60° and 
90° the MS criterion can suitably predict the failure. However, for fiber orientation angles 
between 30° and 60°, a linear combination of the two failure criteria may be adopted as 
described by Eq. 4[20]. 


UF) =GUF ye tQoUP ry (4) 


mixed 


Where (UF) mixeq denotes the mixed index factor. The index factor is the ratio of the stresses to 
strength in which the values of greater than 1 represent failure. a, and a2 are the linear 
combination constants, which can be typically set to 0.9 and 0.1, respectively [20] . UF) ys and 
(IF)rw denote respectively index factors associated with the MS and TW failure criteria defined 
as: 
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oO oO oO o,\o 
IF ),,< = max ! or — |,| ~or— |, 5 
=a (For) oH) : 


UF )ry = Fo,+ Fo, + Fo; Ht E05 a FgOin (6) 


F= ° ae er 
ae oma). Camas aml A 
1 1 
F = 3. aes 7 
© X|X| YK @ 
1 
Fos = 


where X;, X., Y; and Y, denote the ultimate tensile and compressive strengths in the fiber and the 
transverse directions, respectively. S denotes the in-plane shear strength of the lamina. 07, 02, and 
072 are two lamina normal and shear stresses in material coordinates. 


2.2. Composite laminate design for manufacturing guidelines 


Several guidelines have been established based on previous experience from design, test, and 
manufacturing to support the designer utilizing best the composite material strength. The purpose 
of these guidelines is to avoid the occurrence of critical failure modes and simultaneously 
guarantee a manufacturable design. These guidelines are usually defined on the ply level in 
thickness or the longitudinal direction. Bailie et al. (1997) [34] and Kin et al. (2005) [35] 
suggested several rules for stacking sequence design. Some of these design rules are described 
below: 


(1) Symmetry: In most of the real applications, the stack of a laminate are selected to be 
symmetric about its mid-plane to make the design analysis simple. Symmetry also offers 
simplicity in testing, the description of allowable properties, and manufacturing. In symmetric 
laminates components of membrane/bending matrix, as defined in Eq. (3b), are zero. Vannucci 
and Verchery (2001)[36] discussed a more general state for uncoupling of membrane/bending 
deformations in laminated composites. 


(2) Balancing: A balanced layup represents uncoupling in in-plane normal and shear responses. 
Moreover, balancing condition decreases the values of bending-twisting coupling terms in D 
matrix resulted in a simpler bending response of the laminate (see Eq. 3). Montemurro (2015) 
[37] discussed the conditions of a fully orthotropic laminate and uncoupled membrane and 
bending responses. 


(3) 10% rule: A practical laminate should have a minimum 10% of layers should have fiber 
orientations 0°, +45°, and 90° [34] . This choice makes the laminate fiber dominated, and 
provide good damage tolerance and durability [38]. 
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(4) Contiguity constraint: The thickness of a lamina is directly related to the interlaminar stress, 
which affects the occurrence of free edge delamination. This design rule concerns a limitation to 
the thickness of adjacent plies with the same fiber orientation angle. According to Niu (2010) 
[39], for a laminate with a standard ply thickness (0.127 mm), not more than four contiguous 
plies should be used. 


(5) Separation of +@ layers: Separation of layers with +6 fiber orientations can lessen 
interlaminar shear stresses and therefore gives more delamination strength to the laminate [35] . 


(6) Disorientation constraint: when the relative fiber orientation angles difference between two 
successive layer is greater than a limit value, significant interlaminar stresses will be developed, 
which can promote failure. To postpone matrix cracking and the occurrence of delamination, it is 
better to use 45° as a limit for the ply angle difference between two adjacent plies [40]. 


3. Optimization method 


Particle swarm optimization (PSO) is a population-based random search technique motivated by 
the intelligent social behavior of animals. [41]. This method is first proposed by Kennedy and 
Eberhart (1995) [42]. The main idea behind the PSO algorithm is using a swarm to 
simultaneously search a large region in the feasible space coupled with artificial life 
characteristics. 


Consider a set of particles that are distributed uniformly across an n-dimensional search space. In 
each PSO iteration, the position of the particle k is updated by the addition of the vector v;, called 
velocity, to the previous position vector (Eq. 8). The velocity vector is calculated according to 
Eq. (9), in which two vectors representing the difference between the current and the best 
position of the particle, and the current and the global best position of the particle are 
incorporated. 


Kit = Xi +yin (8) 
vil = wv +c,R,(P? -X))+c,R,(P% -X}) (9) 


where X, =[X/, X5, ---, X']', P? =[P’, P’, -, PJ’, and P® =[P’, PB’, ---, P’]’ denote the 


n n 


state of particle k at iteration i, the best state of particle k in the search space, and the global best 
position reached by all particles, P* =[P*, P’, ---, P’]’, respectively. R; and R> are two 
random numbers that vary between 0 and 1, and c; and c2 are cognitive and social parameters 
which are problem-dependent constants. In fact, c; and cz control the magnitudes of particle 
velocity, and hence its movement, along the directions of a particle best position and the global 
best position of the swarm. A limit velocity, Vax, is usually defined to prevent particles from 
moving off the search space. Accordingly, the velocities are bounded as: 


i+] 

v v, 2V 

1 2 2 

v an max a max (10) 
—v vy. Sy 


max k max 


EF. Javidrad et al./ Journal of Soft Computing in Civil Engineering 4-1 (2020) 01-28 9 


@ is a multiplier, known as the inertia factor, used to further control the velocity by making a 
balance between exploration and exploitation. Exploration and exploitation are referred to as the 
ability of a particle movement across the entire search space and inside a local region, 
respectively. The balance between exploration and exploitation is important to approach the 
global optimum point by facilitating jumping out of the local optimum points. Although there are 
many suggestions for this function, a linear step function defined as Eq. 11 is usually used [43]. 


0 gg —{ San =Pa (11) 


where ia, is the maximum allowable solution iteration number. @max and Wmin are two limit 
values of . 


The main shortcoming of the PSO algorithm is that it may trap into a local optimum if an 
improper balancing between exploration and exploitation is provided. Javidrad and Nazari 
(2017) [22] suggested a hybrid method to improve the social performance of a swarm by a 
systematic updating the global best position of the swarm. In this respect, when no further 
improvement in the global best position recognized in a PSO cycle, the SA algorithm is 
implemented to improve the global best position of the swarm. 


SA is a global optimization method based on the Metropolis Monte Carlo procedure. SA uses the 
Metropolis method to accept a state probabilistically with higher function value. The acceptance 
probability is defined as 


(12) 


1 F (Xi) S A(X) 
exp(—Af/T) Otherwise 


Where Afis the variation of the function between two successive points, and 7’ is the acceptance 
parameter called temperature. Temperature controls the exploitation capability of the SA 
algorithm by limiting the acceptance probability. In fact, by sequentially decreasing 7} the 
exploitation ability of the SA is promoted [44]. There are many temperature change functions, 
such as linear, logarithmic, and exponential functions [45]. The most widely used temperature 
decrement rule is: 


T= RT, (13) 
where Rr is a positive decrement constant between 0.8 and 1. 


Referring to the literature, there are a few algorithms developed based on the combination of 
PSO and SA. Xia and Wu (2006)[46] introduced a hybrid PSO-SA method in which PSO 
generates an initial candidate point for the SA analysis. In another approach (Wang & Li, 
2004)[47], the position of each particle within the PSO optimization process is improved by the 
SA process. The integration of the Metropolis acceptance criterion with the PSO is another 
hybrid PSO-SA method, which has been suggested by Shieh et al. (2011) [48]. A more 
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comprehensive literature review on the subject of the hybrid PSO-SA optimization method has 
been given in Javidrad and Nazari (2017) [22]. 


The proposed hybrid algorithm essentially integrates good local search capability of the SA 
optimization method within the PSO. Fig. 2 shows the overall concept of this optimization 
method. 


Generation of particles 


PSO algorithm 


| 


0ud3I9AU0Z, 


Shrink local > 
search region Improvements detected for a 
(decrease limited number of iterations 
temperature) = ] 
jam 


SA local search for P® | 


I 


Update global best point 


Fig. 2. Conceptual model of the hybrid optimization method. 


This algorithm has been implemented and applied to unconstrained single objective optimization 
of several mathematical functions with various dimensions by Javidrad and Nazari (2017) [22]. It 
is reported that including SA within the PSO, upgrade stability and convergence rate of the 
solution process. 


In the present study, however, a modified SA is implemented within the PSO-SA hybrid method. 
In this modified version of SA, the temperature change is adaptively related to the function 
decrease which cause the inner loop of the SA is ignored [49]. It is noted that the inner loop of 
the SA algorithm is related to the length of the Markov chain representing temperature stability 
within each SA loop. 


4. Optimization problem description 


The optimization problem concerns the design of a minimum weight laminate under the stiffness 
and strength constraints and to satisfy conventional design rules. 


The optimization problem is formulated as a bi-objective problem as follows: 


min (W, f) (14) 


where 
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NoL 


W = pAY It, >) 
i=l 

f =max(0,g,(X,)) f=1,2,3 i=1,2,--,n (16) 
MP Sp< — Sp" 7 

g(X,)=>)| ~2—* | -e<0 §=1,2,---,n (17) 
m=1 SP, 

py 20, Fai en (18) 

UF) 1 

MC 

8(X,) = ih (X,) <0, §=1,2,---,n (19) 
k=1 


where n and NoL denote the number of laminae and number of the local regions, respectively. 
, t, A, Sp*, Sp’, MP and, « are the density, lamina thickness, panel area, the required and 
calculated stiffness parameters, the number of stiffness parameters to be included in the analysis 
and the stiffness tolerance value, respectively. (Finixeq); denotes the index factor determined for 
the ith lamina, and (/F),i is the allowable index factor. Similar to the factor of safety, the choice 
of the allowable index factor is dependent on the design and specific application. To provide an 
adequate margin of safety, in this work, the value of 0.83 has been chosen for the allowable 
index factor. This value gives a factor of safety of 1.2 to the failure of layers. MC is the number 
of design rules to be included in the design process, and_X; is the vector of design variables. It is 
noted that for a single lamina, three variables (¢, @, /) are incorporated into the optimization 
formulation. The first two are lamina thickness and orientation angle, which are continuous or 
discrete variables, and the third one is an integer contiguity variable that determines how many 
panels a lamina extended. The contiguity variable can be greater than one only for problems 
concerming blended laminate layup. 


X,=(t,0,),,  i=l,---,n 
-90 <A<90° 

t,<t<t, 

ISIS Wipes 


(20) 


Where t, tu, and Neg are the lower and upper limits of lamina thickness, and the number of 
regions in the laminate, respectively. i, is an integer function introduced to characterize the 
degree of the manufacturing constraints violation. 4, varies between 0 for all constraints violated, 
and 5 for all constraints satisfied. The symmetry condition is always enforced in this analysis and 
thus is not included in Ay. 


(21) 


O if Athconstraint is satisfied 
« "1 if kth constraint is not satisfied 
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5. Implementation 


e Initial swarm is generated and distributed uniformly across the design space (i.e., within the 
lower and upper bound of the design variables). For all case studies discussed in this article, the 
population size was 50. 


e The parameters of the optimization analysis have been selected as in Table 1. It has been 
shown that these are optimized parameters for a range of mathematical and engineering problems 
[22]. These parameters remained constant in all case studies presented herein. 


Table 1 

The Optimized parameters of the proposed hybrid PSO-SA method. 
To Rr max Omin Ci C2 
1.0 0.95 0.55 0.1 0.8 1.7 


e For each particle, the best position (P”) at iteration i are determined. Similarly, the global best 
position (P*) of the swarm is calculated. In the determination of these two positions, the 
following selection criteria are implemented: 


(1) If the value of fis zero for two particles, the particle with the lower W value is selected. 


(2) If the value of fis zero for one particle and nonzero for another, the particle with zero fis 
selected. 


(3) If the value of fis not zero for two particles, the particle with the lower f value is selected. 


At each iteration, every P® with zero value of violence function f can be a potential optimum 
solution, i.e., Pareto solution. These solutions are stored for the final selection process. 


e The velocity and position of each particle are updated according to Eq. (8) and Eq. (9). The 
new swarm is designated by X*"’. 


e The procedure continues until no improvement found in the global best solution during a PSO 
cycle. Then the method shifts to the SA, and the solution procedure continues for P®. The SA 
loop length is pre-defined, and the final P* determined during this loop returns to the PSO. The 
SA method used here is a modified one in which temperature is only updated when a move is 
accepted. By using this modified SA method, only one loop is necessary. 


e The following convergence criteria have been used to terminate the solution procedure: 
(1) The number of iterations reaches the prescribed maximum value. 
(2) The temperature reaches to a tolerance value (namely, 10°). 


e When the solution process terminates, the point with minimum VW is selected from the solution 
vector as the final optimum solution. 
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6. Numerical results 


6.1. Design of a laminate with maximum 4A), stiffness coefficient 


In this example, we seek ply stacking sequence of a laminate with n continuous plies that 
maximize A». coefficient for a graphite-epoxy laminate. The composite is assumed to be 
symmetric and balanced. This problem, which is a variation of a composite design problem 
undertaken by Grosset et al. (2005) [50], has been approached by Venter and Haftka (2010) [51]. 
A constraint on the effective Poisson’s ratio (Ver) and limitations in fiber orientation angles and 
ply thicknesses have been considered. The optimization problem can be stated as 


Maximize Ag) =h(U, UV, +U3V3) 
Subjected to 0.48 < Ve" < 0.52 
a5 A.< 5° or 
40 <@<50 or (22) 


85 <0, <95° 
0.001< t, < 0.05 


where 
> h/2 2 n 
V3) = re | {cos 20, cos 40} dz = — Dd talcos26,, cos 46, } (23) 
0 k=l 
A U4 —U3V3 


14 = ~ * * 24 
a Ayo U; _ UV, + U3V3 ( ) 


and G and t, are ply angles and thicknesses (in inches), respectively. The U; values are material 
invariants as given in Table 2. 


Table 2 

Invariant parameters of the graphite-epoxy composite. 
Parameter Value (psi) 
U; 0.8897x107 
U; 1.025410! 
U; 0.2742x107 
Us 0.3103x107 


Choosing n=3, we have only six design variables. This problem has been examined by a variety 
of methods such as single objective PSO algorithm using exterior quadratic penalty function and 
the adaptive penalty method, multi-objective PSO (MO-PSO), The modified multi-objective 
PSO using the crowding distance as the second selection criterion, and the PSO-SA algorithm 
presented herein. A maximum of 3000 function evaluations and a temperature limit of 10° were 
set as the termination criteria. To attain statistical results, several independent optimization runs 
were performed. The results are summarized in Table 3. The success rates, which is defined as 
the ratio of successful runs to the failed runs, are also given in Table 3. The successful runs are 
those reached to the known global solution. 
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The calculated layup stacking sequence by the proposed bi-objective PSO-SA method is very 
close to that of the other methods. When compared to the other multi-objective formulations as 
given in Table 3, the calculated maximum A.) is close enough to the global solution, 1.25*10° psi 
[51]. Typical variations of the A». and ve are shown in Fig. 3 and Fig. 4, respectively showing 
that the method explores the design space thoroughly and then stably approaches to the global 
optimum point. 


Table 3 
Optimization results for 4. maximization problem. 
Single objective PSO" . Modified PSO-SA 
Parameter Quadratic penalty Adaptive penalty MO-PSO one (present wok) 
Best layup 
orientation [+95, +44.3, 444.5],  [+95, 444.3, 444.5], [+95, +43.7, +42.3], [+95, +44.8, 444.4],  [+95,4+43.3, +44.4], 
(degrees) 
Thickness (in) [0.0304, 0.05, 0.05], — [0.0304, 0.05, 0.05], (0.0323, 0.05, 0.05], [0.030, 0.05, 0.05], [0.0312,0.05,0.05], 
Vert 0.4800 0.4800 0.4800 0.4800 0.4805 
Success rate (%) 97 37 29 100 100 
Maximum A) 1.2505 x 10° 1.2506 x 10° 1.2434 x 10° 1.2503 x 10° 1.2485x 10° 


* Taken from Venter and Haftka (2010) [51]. 
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Fig. 3. Convergence behavior of Az. by the hybrid PSO-SA optimization method. 
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Fig 4. Convergence behavior of Veg by the hybrid PSO-SA optimization method. 
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As a continuation of the problem, we implemented discrete fiber angle orientation, layer 
thickness, and design rules as described in Section 2.2 to find a better design. It is assumed that 
the laminate comprises 32 layers with the layer thickness of 0.0083. This assumption gives a 
weight near to the previous results. The results for two fiber orientation angle steps (5° and 15°) 
are summarized in Table 4. The results given in Table 2 show that including more design rules to 
the optimization analysis will affect the maximum value determined for A22. By reducing the 
design rules constraints, the maximum value of Ay) can be greater than 1.25x10°. 


Table 4 
Optimization results for Az. maximization problem including discrete variables and design rules 
constraints. 


Angle Case" Layup Ax x10° Veft 
1 [55,-50,-55,-50,45,-55,-65,-55,50,55,50,-45,55,652]s 1.385 0.4802 
5° 2 [65,-60,-50,55,45,55,-50,-60,-65,60,50,-55,-45,-55,50,60]s 1.385 0.4822 
3 [-60,-55,-50,-605,-30, 15,60,55,50,602,30,-15,-60]s 1.210 0.4809 
4 [-55>,-50,-40,0,45,50,90,555,50,40,0,-45,-50,-90]s 1.177 0.4809 
1 [-602,45,-45,-605,60,-45,60>,-45,45,602,-60,45]s 1.359 0.4886 
A 2 [-60,-45>,60,-60,-45,-60,60,455,-60,60,45,60.]5 1.359 0.4886 
3 [-45,-60,-45,-603,-45,0,45,60,45,603,45,0]s 1.171 0.5084 
4 [-60,-453,0,45,90,45,60,453,0,-45,-90,-45] 1.109 0.5165 


“1: Including all design rules except the 10%, Contiguity and Disorientation design rules. 
2: Including all design rules except the 10% and Disorientation rules. 

3: Including all design rules except the 10% rule. 

4: Including all six design rules. 


6.2. Layup sequence design of a laminate under in-plane forces and moments 


As the second design case study, a laminate made of Kevlar 49/epoxy with two different load 
cases has been considered. The laminate is assumed to be consist of twenty-ply groups. The 
optimization problem is to find the minimum laminate thickness based on design rules, to have 
the prescribed effective engineering coefficients and not to fail under the given set of loads. 


To compare the results of this work with the published data, we set the material properties, load 
cases, and the prescribed engineering coefficients the same as those given by Javidrad et al. 
(2018a) [51] (see Tables 5-7). The Fiber orientation angle of each lamina was assumed to be 
multiples of 5°, 15°, and 45°. Two termination criteria were used: (7) the number of function 
evaluations exceeding 1x10‘ and (ii) the temperature approaching less than 1x10°. 


Table 5 
Typical elastic material constants of Kevlar 49/epoxy composite. 
Ey Ey Gi X, X, Y, Y, S 


Vi2 


p 
(GPa) (GPa) (GPa) (MPa) (MPa) (MPa) (MPa) (MPa) (kg/m) 
76.0 5.5 2.3 0.34 1400 235 12 53 34 1450 
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Table 6 
Loading conditions. 


Loading condition (1) Loading condition (2) 


N, (N/m) 20000 2000 
N, (N/m) 130000 1400 
My (N/m) 0 0 
M,. (N.m/m) 0 0 
M, (N.m/m) 0 150 
M,, (N.m/m) 0 0 
Table 7 
The required effective engineering constants. 
Ee Ey Gry yo EX, Ey, Cs Pe 
(GPa) _(GPa)___(GPa) (GPa) (GPa)___—(GPa) a 
26.0 19 14 0.55 16 20 10 0.35 


Superscript 0 and f stand for in-plane and flexural properties, respectively. 


The best layup sequences found for each set of loads together with the results of the single- 
objective optimization process of Javidrad et al. (2018a) [51] are summarized in Table 8. The 
best layup sequences are those with zero function f (Eq. 16) and the least weight, which is 
selected from all potential optimum solutions. A comparison with this reference indicates better 
or similar results for all cases except case no. 6. In case no. 6, however, our index factor is lower 
than that given in the reference data. The effective engineering constants, as given in Table 9, 


show a minor difference with the required effective properties. 


Table 8 
The best calculated layup sequences. 


Fiber orientation angle 


Fone Case Stacking sequence step Wee Ieee 
condition (kg/m*) factor 
(Degree) 
1 [45,90,45,403,0,-45,-90,-45,-403,0], 5 5.278 0.824 
als [45,90,45>,353,0,-45,-90,-45,- 5 6.032 0.683 
353,01. 
2 [-45>,-90,-45,- 15 5.278 0.822 
(1) 302,0,452,90,45,302,0], 
a [-45,-90,-454,- 15 6.786 0.563 
302,0,45,90,454,302,0], 
3 [-45>,-90,-453,02,452,90,453,02], 45 6.032 0.700 
3° [-45>,-90,-453,02,452,90,453,02], 45 6.032 0.700 
4 [-453,-90,-50,- 5 6.032 0.741 
255,0,453,90,50,25>,0], 
4° [45,90,454,25>,0,-45,-90,-454,- 5 6.786 0.690 
252,01; 
(2) 5 [-453,-90,-455,- 15 6.032 0.828 
15,0,453,90,45>, 15,0], 
5" [45>,90,454,02,-45>,-90,-45 4,02], 15 6.786 0.633 
6 [-45>,-90,-454,02,452,90,454,02], 45 6.786 0.622 
6 [452,90,453,0,-452,-90,-453,0], 45 5.278 0.760 


* represents the best results obtained from Javidrad et al. (2018a) [20]. 
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Table 9 
Effective engineering constants determined for each layup. 
E° E° Go , E! EL Gf / 
Case i é : vy, ~ i v., 
(GPa) (GPa) (GPa) - (GPa) (GPa) (GPa) 
1 24.130 20.090 14.291 0.533 15.476 22.634 9.960 0.306 
2 28.248 20.062 13.291 0.523 17.267 20.823 10.078 0.338 
3 28.470 21.181 12.984 0.491 18.440 20.562 10.653 0.351 
4 28.145 20.178 13.292 0.520 17.251 19.118 10.347 0.369 
5) 26.814 20.669 13.520 0.511 16.719 18.805 10.642 0.382 
6 26.473 20.232 = 13.698 0.524 17.239 19.687 10.719 0.370 


6.3. Design of plate under bending loads and blending constraints 


The optimum design of a composite plate with layup blending constraint is considered as the 
final example. In general, multi-region non-blended panel weight optimization with no 
continuity of laminate layup may be structurally non-profitable due to discontinuities. Although 
imposing the layup blending constraints gives rise to the complexity of the optimization process, 
it yields a more practical design. There are many approaches addressed in the literature for the 
design of composite laminates using layup blending concepts such as greater—than-or-equal-to 
blending [52] , sheared-layered blending [53], guide-based blending [54], and stacking sequence 
table based blending [55]. 


In the greater-than-or-equal-to blending method, a key panel is defined from which each ply 
emanates. The key panel is often the most loaded or the thickest region. Although each ply or 
lamina is allowed to be discontinued along the panel borders but once dropped, it is not allowed 
to be combined back into the other panels. Although this method is quite simple and can easily 
be implemented to the design process, it is not suitable for designing structures with several key 
panels. 


To date, most of the works in composite structures design by layup blending methods have been 
centered on the single objective optimization using GA. In this example, however, we aimed at 
designing a cantilever composite plate using the concepts of the greater-than-or-equal-to 
blending by the proposed PSO-SA optimization in the framework of the bi-objective 
formulation. The plate is shown in Fig. 5, which is decomposed into three local panels. 


As stated in Section 4, for a single lamina, three variables (¢, @ /) are incorporated into the 
optimization formulation. / is an integer contiguity variable that determines the extension of the 
lamina through the plate. In this example, we divided the plate into three equal panels, in which 
the Panel | is assumed to be the key region. In this respect, /=0 denotes that the lamina does not 
exist, /=1 designates the lamina only exists in the Panel 1, /=2 denote the lamina extended 
through the Panels | and 2, and finally, /=3 represents the lamina continues through all panels. 


The design optimization problem is stated as a minimum weight design of the plate subjected to 
the maximum global deflection (in the Z-direction), first ply failure, local ply contiguity, and 
conventional design rules. 
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Fig. 5. Geometry of the cantilever composite plate. 
In this example, the stiffness of the plate is represented by the global tip deflection in the Z- 
direction. In this regard, Eq. 17 is substituted by 


~1/<0.01 (25) 


| | Wrax 


Ww 


where w,,, and w denote the maximum calculated deflection in the z-direction and the 


maximum allowable deflection, respectively. 


The in-plane material properties of the plate, which is made of carbon fiber epoxy laminate, are 
given in Table 10. The through-thickness material properties are assumed to be identical to those 
of the in-plane properties. A uniform distributed load of 0.1 kN/m is applied to the tip of the 
plate, as shown in Fig. 5. 


Table 10 

Typical elastic material properties of carbon/epoxy laminates. 
Ey Ey Go, X Xe Y Y. S p 
(GPa) (GPa) (GPa) ! (MPa) (MPa) (MPa) (MPa) (MPa) kg/m’) 
168 9.5 46 0.27 826 118 37 27 20 1600 


A finite element uniform grid (25x15) consisting of 312 4-noded S4R composite shell element of 
ABAQUS software has been used to model the plate behavior. This shell element is well suited 
for moderate to thick composite plates since it accounts for transverse shear flexibility. The stress 
distribution through the thickness of each layer is determined by specifying several integration 
points through the thickness. 
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The optimization process has been linked to the finite element software to determine the 
maximum deflection and stresses as needed for the design optimization procedure. All 
parameters for the PSO-SA optimization analysis were set as given in Table 1. The maximum 
iteration was 9000. The results for 0.5 < JF < 0.83 and w=100mm are summarized in Table 11. 


Table 11 
The optimum blending layup sequences for all three panels. 
Solution : Index Weight Maximum Deflection 
No. panel Stacking sequence Factor (kg) (ni) 
1 [-455,-15,-30,03, 45>, 15,30,03], 0.592 
1 2 [-45,-15,03,45, 152,03], 0.564 7.6 0.1073 
3 [-45,03 452,02], 0.513 
1 [-455,-15,-25,03,453,15,25,03], 0.591 
2 2 [-455,-152,03,452 15,03], 0.553 7.6 0.106 
3 [-453,03,452,03], 0.506 
1 [-45>,-15 ,-30,03,45>, 152,30,03], 0.592 
3 2 [-455,-15,-30,03,45,15,30,03], 0.560 8.0 0.0986 
d [-452,03,452,03]. 0.511 
1 [-45>,-152,-45,03,459,15>, 45,03], 0.595 
4 2 [-455,-15,-45,03,459 ,152,45,03], 0.576 8.0 0.1009 
7 [-452,03,452,05]. 0.519 
1 [-453,-20,-53,0,453,20,53,0], 0.514 
5 2 [-453,-20,-53,0,453,20,53,0], 0.564 8.4 0.101 
2 [-453,-20,0,453,20,0], 0.518 
1 -45,-50,-15>,- 
sino en mee! 
° z [-45,-152,-204,0,45,152,20,,0]. 0.519 oe O.087 
3 [-45,-155,0,45,155,0], 0.512 
1 [-455,-152,-454,0,459,15,454,0], 0.581 
7 2 [-455,-152,-454,0,459,152,454,0], 0.666 8.8 0.0986 
2 [-45,-155,0,45, 15,0], 0.600 
1 [-453,-15,-354,0,453, 15,354,0], 0.510 
8 2 [-453,-15,-354,0,453, 15,354,0], 0.737 9.2 0.107 
é [-453,-15,0,453,15,0], 0.645 


Among the results given in Table 11, the best solution is the solution designated by No. 1, which 
can be regarded as the global optimum solution. It is noted that if the maximum number of 
iterations increases, better solutions may be found. All design rules except the 10% rule have 
been included in the analysis. When the 10% rule applied to the analysis, no feasible solution 
was determined. The variations of two objective functions, together with the optimum solution 
points are plotted in Fig. 6 in which the solutions with zero fare highlighted by solid bullets. It is 
seen that after finding a solution, the process explores the search space further to find other 
solutions. 
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Fig. 6. Variation of objective functions with respect to the solution iteration. (Solid bullets indicates the 
optimum points) 


7. Conclusions 


An optimization procedure based on a hybrid PSO-SA method has been presented for the design 
of laminated composites. It has previously been shown that the hybrid PSO-SA optimization 
method provides good exploration and exploitation behavior, which undoubtedly is significant 
for solving constrained problems with narrow search spaces. Stiffness, strength, continuous 
layup blending, and several design rules have been included in the analysis to enhance the 
applicability of the results for actual engineering problems. This formulation combines all 
constraints into an objective function and therefore has less complexity as compared to general 
single-objective constrained optimization approaches. Moreover, several optimum solutions may 
be determined for a particular problem, which increases the flexibility of the design. Three 
examples, including a composite plate bending problem with layup blending constraint, have 
been presented to show the effectiveness of the proposed method. It can be concluded that the 
proposed method describes an effective and relatively simple process for the design of laminates, 
and it can be used as a basis for the development of an sophisticated system for design of 
laminates. 
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8. Replication of results 


In this section, we provide a pseudo-code of the PSO-SA optimization method and some 
benchmarks for the validity of the proposed method. As stated previously in this article, the PSO- 
SA method is a hybrid stochastic approach that integrates the good local search of the SA within 
the PSO. In this regard, improvements in the global best solution will be performed by the SA 
method, if certain conditions were met. A new version of the SA is implemented here in which 
temperature reduction is adaptive [49]. In the present pseudo-code, the parameters must be set 
according to those given in Section 5. In Section 5, we also explained the implementation of the 
proposed method in solving problems in composite materials design optimization. 


Procedure PSO-SA 
n @ NumberParticles //{set initial parameters} 
max “ MaxLimit 
min “ MinLimit 
imax *” MaxIterationNumber 
c, — CognitiveParameter 
C) + SocialParameter 
Vmax * UpperVelocityLimit 
Nmax * SALoopLength 
. 79 — SAInitialTemperature 
. Ry  SATemperatureDecrementParameter 
. // {generate initial state} 
. for i=1 ton // {initialize n particles} 
X* — GeneratelnitialParticles 
Ph, — X* 
. Endfor 
. T=T 
. // {set Peas the best of the initial population} 
. Stopflag=.false. 
. // {determine point of best solution} 
fer”) & min (f(X) 
pee Xhmin 
. fori from | to ina, //{main loop} 


SO DOSS ON SIN SY OR 


Se eee 
wWNrF © 


Se etl tll loo 


NNNNW 
RwWNr OS 


®  Dmax-( @max- Omin) *t/ imax 
for k from | ton 
R, © Random(0,1) 
R, + Random(0,1) 
. // {update velocity and position of each particle} 
ve & ov! t ex Ry (Phy -X4)+ Co Ry (Pg'-X;) 
xi + Xi+y,it! 
. //{check velocity for upper limit} 
If (v aaa GEWVmax) V Sie 
If (v ge oe ee | Vg Si 
. // {update best of each particle} 
P bist" =P bi 
. // {compute the global best solution of the swarm, Pg't 
. // {if Pg' is not improved in a PSO cycle, then send it 
. // to the modified SA procedure for an improvement} 


WWWWWWWWWNNNDND LY 
SAIAAARKRWNHESOANIAAUA 
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39. // {Modified SA Procedure} 


40. XSA;+ Pg’ 

41. for j from | to Max //{SA loop) 

42. XSA, + PickRandomState 

43. // {compare objective function values} 

44, if (f(XSA,)=< f(XSA,) then 

45. XSA; — XSA, //{accept move} 

46. T — Tx Ry //{temperature update} 
47. Else 

48. //{check for Metropolis criteria} 

49. R + Random(-1,1) 

50. if (R<exp(-f(XSA,)+ f(XSA,))/T) then 
51. XSA; + XSA, //{accept move} 
52. T + Tx Ry //{temperature update} 
53. Endif 

54. Endif 


55. //{update global best} 

56. Pg' — XSA; 

37. endfor 

58. return // {to PSO} 

59. // {end SA procedure} 

60. //{check for stop criteria} 

61. If Stopflag=.true. then return 
62. endfor 

63. end 


To verify the developed computer code and test the performance of the proposed bi-objective 
formulation, four mathematical constrained optimization problems have been resolved using the 
proposed bi-objective optimization analysis. The test functions and constraints are given in Table 
12. Kim and Myung (1997) [56] described seven algorithms for mathematical constrained 
optimization problems based on two evolutionary programming (EP) methods. For comparison 
purposes, we selected two of these algorithms for each test function, which resulted in the best 
and the worst solutions. The obtained optimum solutions, together with the results from the 
literature, are summarized in Table 13. It is noticeable that the proposed bi-objective algorithm 
presents good results when compared with these reference values. The convergence curves are 
shown in Fig. 7 through Fig. 10. 


Table 12 
Test functions. 
Function Objective Constraints Bounds 
2 
x, +25 20 0.5 < x, < 0.5 
Fl 100(x, — x7)? + —x,)? ei ’ 
xX + x7 20 0<x <1 
iy Oe eBay 8x 22 O<x, <3 
F2 TX — XQ 
Xy — 4x} +3247 —88x7 + 96x, < 36 O<x. <4 
F3 2 XyX_ — 2520 2< x, <50 
mane x? +x -2520 0< x, <50 
Axf 20 
F4 (4, - 2)? + (a -1)? 5 


Xp +X <2 
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Table 13 
Comparison of the results for mathematical constrained optimization problems. 
Function Method’ Best Mean Worst 
Algorithm #4 0.25 0.25 0.25 
Fl Algorithm #7 0.25 0.25 0.25 
Present 10077) 0) 0) 
Algorithm #1 -5.50801 -5.50801 -5.50679 
F2 Algorithm #5 -5.50801 -5.50801 -5.50801 
Present -S.50801 (2. wre ree 
Algorithm #2 5 5.12323 5.49622 
F3 Algorithm #5 5 5.0002 5.00082 
Present See ene 
Algorithm #3 1 1 1 
F4 Algorithm #7 1 1 1 
Present Pee ee 


“The results designated by Algorithm #1 thru Algorithm#5 and Algorithm#7 were taken from Kim 
& Myung (1997) [56]. 
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Fig. 8. Convergence curves for the function F2. 
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